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This paper discusses a computationally efficient algorithm for estimating the safe ma- 
neuvering envelope of damaged aircraft. The algorithm performs a robust reachability 
analysis through an optimal control formulation while making use of time scale separation 
and taking into account uncertainties in the aerodynamic derivatives. This approach differs 
from others since it is physically inspired. This more transparent approach allows inter- 
preting data in each step, and it is assumed that these physical models based upon flight 
dynamics theory will therefore facilitate certification for future real life applications. 


N omenclat ur e 


a angle of attack [°] 

f3 angle of sideslip [°] 

7 flight path angle [°] 

A uncertainty [-] 

A uncertainty vector [-] 

p air density [kg/m 3 ] 

(f> state trajectory 

ip bank angle [°] 

X course [°] 

C m aerodynamic derivative [-] 

D drag [N] 

f nonlinear continuous function 

Fa. aerodynamic forces in the aerodynamics frame of reference [N] 
g gravitational acceleration [m/s 2 ] 

1 Invariant set 

K set of states 

L lift [N] 

L set of states 

l implicit function 

m mass [kg] 

p costate vector 
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Pi 1 P 2 costates 

q dynamic pressure [kg /ms 2 ] 

1 Z Reachability set 

S wing surface area [m 2 ] 

standard deviation [-] 

T thrust [N] 

T time horizon [s] 

t time [s] 

U the set of Lebesgue and bounded measurable input functions 
u input vector 

V Viability set 

V airspeed [m/s] 

V± value function for viability 

V 2 value function for invariance 

W weight [N] 

X\]X 2 state 

x state vector 

x 0 initial condition for state vector 

Yaero aerodynamic sideforce [N] 


I. Introduction 

S afety is of paramount importance in all transportation systems, but especially in civil aviation. There- 
fore, in civil aviation, many developments focus on improving safety levels and reducing the risks that 
critical failures occur. In a recent study by CAST/ICAO, it can be observed that loss of control in flight 
(LOC-I) is the most frequent primary accident cause. This study is based on a statistical analysis of air- 
craft accidents between 2002 and 2011, and indicates that this category accounts for as much as 23% of all 
fatal aircraft accidents and involves most fatalities? Benefit can be gained by developing technology which 
prevents these LOC-I accidents. From a flight dynamics point of view, with the technology and computing 
power available on this moment, it might have been possible to recover some of the aircraft in the accident 
category described above on the condition that non-conventional control strategies would have been ap- 
plied. These non-conventional control strategies involve the so-called concept of fault tolerant flight control 
(FTFC), where the control system is capable to detect and adapt for changes in the aircraft behavior. 

One FTFC strategy option is using a model based control routine. Previous research focused on a physical 
model approach. 2 In this setup, experiments have shown that not only a reconfiguring controller is needed, 
but also some form of flight envelope protection, which prevents the airplane from leaving the safe flight 
envelope and loosing control in flight P The main challenge in this context is determining the new bounds 
of the safe flight envelope after failure, which are then used by the envelope protection algorithm. 3 

Alternative and complementary research approaches for the purpose of loss of control prevention and 
prediction are among others passive adaptive control,® data-based predictive control 5 and real-time optimal 
envelope limit estimation. 6 

Determination of the flight envelope has been done in the literature through various methods. The most 
straightforward methods include wind tunnel testing, flight test experiments and high-fidelity model-based 
computation of attainable equilibrium sets or achievable trim points,®®^ possibly with bifurcation anal- 
ysis® 12 or a vortex lattice algorithm combined with an extended Kalman filter! 13 More complex methods 
include formulating flight envelope estimation as a reachability problem and solving this with level set meth- 
ods and Hamilton- Jacobi equations)® 15 16 17 18 possibly with time scale separation 19 or semi-Lagrangian 
level sets!® Alternative methods rely on linearization and region of attraction analysis)- 21 determining con- 
trollability /maneuverability limits in a quaternion-based control architecture 22 or robustness analysis for 
determination of reliable flight regimes!® An approach suggested by Boeing uses Control-Centric Modeling, 
dynamic flexible structure and load models! 24 In the frequency domain, stability margins can be estimated 
in real time via nonparametric system identification. 25 A complete different alternative is relying on the 
bio-immune system metaphor in combination with artificial intelligence J 26 32l][m More focused techniques 
inspired by flight dynamics exist as well, such as determining the minimum lateral control speed P 9 Online 
methods are preferable from an operational point of view. For fault tolerant systems, the adaptive flight 
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envelope estimation must take into account damage information which will possibly restrict the safe envelope 
significantly. 

From the perspective of the physical approach, the preferred interpretation of the safe maneuvering enve- 
lope considers reachability from the trim envelope. The stable and controllable trim envelope is considered 
an a-priori safe set. The backwards reachable set is defined as the set of states from where (at least one point 
in) the trim envelope can be reached. The forwards reachable set is defined as the set of states which can 
be reached from (at least one point in) the trim envelope. Then the safe maneuvering flight envelope is the 
cross section between the forwards and backwards reachable sets. This interpretation is illustrated in Fig. [lj 
In addition to the safe envelope, the backwards reachable set is considered as the survivable flight envelope. 
After an upset due to damage, turbulence, a wake encounter etc., it is possible to bring the aircraft back to 
a safe trim condition as long as the current flight condition is situated inside the backwards reachable set. 



a-pnon sale & 

= trim envelope 


safe operating set = safe flight envelope 
backwards reachable set \ 

i = sun ' ivable f,i 3jiL enve '°P e \ forwards reachable set 


Figure 1. Safe maneuvering envelope as intersection between forwards and backwards reachability, source: 
van OortBID 


The aim is to perform a combined forward and backward reachability analysis from the trim envelope as 
efficiently as possible, for on-line implementations. Based upon previous research, 6 ! level set methods are an 
excellent candidate. Two of the major challenges are the computational load and how to cope with nonlinear 
systems with higher dimensions. In general, an increase in technology readiness level (TRL) is envisaged. 

Nonlinear systems with higher dimensions can be simplified by considering the principle of time scale 
separation. 19 The structure of time scale separation is analogous as applied for the fault tolerant control 
algorithm developed earlier. 2 The overview can be found in Fig. [2j which illustrates that a nine dimensional 
nonlinear problem is decoupled in three consecutive three dimensional optimization problems. 
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Figure 2. Separation of dynamics over high bandwidth, middle range and low bandwidth 
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II. Optimal Control Formulation 


It has been shown in the literature that maneuvering envelope estimation through reachability can be 
reformulated in the optimal control framework. 14 ^ Consider a continuous time control system: 

x = f(x,u) (1) 

with x G M n , uGf/C M m , f (•,•): M n x U -A M n , a function: 

/ (•) : M n -A M (2) 

and an arbitrary time horizon T > 0. Let U\t^\ denote the set of Lebesgue and bounded measurable 

functions from the interval [£,£'] to U. Define (j) (r, £, x, u (•)) as the state trajectory. Given a set of states 

K C M n , a number of questions can be naturally formulated regarding the relation between the set K and 
the state trajectories <f> of Eq. © over the horizon T. Problems of interest include the following: 

Viability: Does there exist a u G l I[ o,t] f° r which the trajectory (j) of the state x satisfies x G K for all 
t G [0, T]? 

Invariance: Do the trajectories 0 of the state x for all u G 0 j t] satisfy x G K for all t G [0, T]? 

Reachability: Does there exist a u G W[o,t] and a t G [0 , T] such that the trajectory (j) of the state x 
satisfies x G if? 


More generally, the following sets can be characterized: 


V(t,K) = {x G r|3uGW M) Vre [i,T] ,0(r,i,x,u(-)) 
Z(t,K) = { x £ R n I Vu G W [f>T ] , Vr e [t, r] , (/) (r, t, x, u (•)) 
K(t,K) = { X e 1"| 3u e 3r £ [f, T] , ^ (r, f, x, u (•)) 

e K } 
e K } 
£ K } 


By comparing the sets, it can be stated that: 



X(t,K) C V(t,K) C 71 (£, K) 


(3) 

The principle of duality results in the relationship, as illustrated in Fig. [3j 



K(t,K) = (l(t,K c )) c 


(4) 



Figure 3. Illustration of the principle of duality. Left 7 Z(t,K), right Z(t,K c ). 


The above set characterizations can be rewritten as optimization problems. First a connection is estab- 
lished between viability and the SUPMIN optimal control problem. Assume that the set K is related to the 
zero level set of a continuous function l : M n — )> M by iC = {xG M n | l (x) > 0}. Then, it can be stated that: 

V (t,K) = < x G M n | Vi (x, t) = sup min l (0 (r, £, x, u (•)))> 0 > (5) 

1 u { .)eu [t , T] re[t,T\ J 
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The connection between invariance and the INFMIN optimal control problem can be established by 
considering a closed set L, that can be written as the level set of a continuous function l : W 1 -A M, i.e. 
L = {xG R n | l (x) > 0}. Then, it can be stated that: 


X(£, L) = < x G 


|r 2 (x,i)= inf min x,u (■)))> 0 

u(-)eU[ tjT ] r£[t,T] J 


(6) 


These optimization problems have been transformed in Hamilton-Jacobi-Bellman Partial Differential 
Equations (HJB PDE)JUl Viability can be written as follows: 


dVi 

~dt 


(x,t) 


min 

re[t,T } 


sup 

1 i')£U[t,T] 


^ (x, t) f (x, u) 


= 0 


(7) 


Furthermore, for backward integration V\ (x, T) = l (x) holds. Forward integration implies V\ (x, t) = l (x). 
Simlarly, invariance is represented by the following HJB PDE: 


dV 2 

dt 


(x, t) + min 

re[t,T ] 


< inf 

lu(-)GW[ t)T ] 


dV 2 

<9x 


(x,t)f (x, u) 


= 0 


(8) 


In addition, V 2 (x, T) = / (x) holds for backward integration, while V 2 (x, t) = l (x) applies in case of 
forward integration. These HJB PDE’s can be solved by level sets, for which a toolbox is available in 
Mat lab® P 


III. Approach 

The approach to calculate the safe maneuvering envelope is based on the following steps: 

• Define reference boundaries for E, 7 and ip as well as grid step size AE, A 7 and Acp 

• Define an implicit function accordingly over E and 7 . This needs to be done for every value of ip in 
case speed and flight path angle are bank angle dependent, i.e. V = V{<p) and 7 = 7 (ip). 

• The Level sets method toolbox 15 relies on the Hamiltonians in Eq. © and ([5]) as a gradient to evolve 
the implicit function and thus reference boundaries over time 

• The cost functions V\ and V 2 become three dimensional functions, where cross sections reflect the 
situation for a specific time instant t a . 

• A dissipation function is needed to guarantee numerical stability during these calculations. As a 
consequence, slightly conservative results are obtained for the boundaries, but analysis has shown that 
this dissipation has a minor effect on the results. In this specific context, the chosen dissipation function 
is a Lax Friedrichs dissipation function . 15 

IV. Application example 

To illustrate how the envelope estimating algorithm works, a nonlinear 3D aircraft example is considered. 
At this point, only the slow dynamics as specified in Fig. [2] are considered. Future work will extend to the 
faster dynamics. The data used in this example are based on the RCAM (Research Civil Aircraft Model) 
simulation modelP^ The acting forces on the aircraft are illustrated in Fig. [4] for a symmetric flight condition. 


For the complete 3D situation, the equations of motion are written as follows P 

Fa x ~ W sin 7 = mV (9) 

F a z cos ip + Fa y sin ip + W cos 7 = —mV 7 (10) 

— Fa z sin ip + Fa y cos ip = mV x cos 7 (11) 

Where the aerodynamic forces can be simplified assuming small aerodynamic angles a and /3: 

F Ax = Tcos (3 cosa — D (V,a) & T — D (H, a) (12) 

Fa z = — T sum — L (E, a) ~ — L (E, a) (13) 

F Ay = —T sin (3 cos a + E aero ( V, /?) ~ E aero (E, /?) (14) 
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n\y* 



Figure 4. 


Acting forces on the aircraft model, source: 


Lygeroi^ 


with the following expansions for lift L, drag D and sideforce Ya ero : 


D(V,a) = qS(C Do +C Da a + C Da2 a 2 ) (15) 

L(V,a) = qS(C Lo +C La a ) “ (16) 

Y aem (V,(3) = qS(C Y 9 0 ) (17) 

where the dynamic pressure q = l/2pV 2 . 

The corresponding numerical values are: m = 120 • 10 3 kg, g = 9.81 m/s 2 , W = mg, p = 1.225 kg/m 3 (sea 
level), S = 260 m 2 , C Lo = 1.0656, C La = 6.0723, C Do = 0.1599, C Da = 0.5035, C Da2 = 2.1175, C Yp = -1.6. 


In the perspective of reachability from stable and controllable trim conditions, the primary states of 
interest are airspeed V and flight path angle 7 . Considering time scale separation as presented in Fig. [ 2 j the 
virtual inputs for the slow dynamics are roll angle p, angle of attack a, sideslip angle (3 and thrust T. This 
framework and combining Eqs. (f9|)- (p!7|) allows to define the general dynamics f in Eq. 0 by the following 
differential equation: 


V 

7 


pS \r2ri 
2rn 


D 0 ~9 sin 7 


- y COS 7 


cos a cos /3 

(cos p sin a cos (3 — sin p sin f3) y 


T 

m 


(< C Da a + C Da2 a 2 ) 

4- 

0 

. £V(C Lo +C La a) cosp _ 


.-2^ VCY ^ Sin ^. 


Assuming small aerodynamic angles, as earlier, simplifies the differential equation: 


V 


~-£v 2 C Do -gsm 7 ' 

A 

v 

m 

T + 

~2§i v2 ( C D a cx + C Da 2 a 2 ) 

A 

0 

i _ 


— ^ cos 7 


0 

. 2^ V ( C Lo + C L a Ot) COS P _ 


r &VC Ye (}8mv 


(18) 


(19) 


This results in a Hamiltonian function with decoupled virtual inputs T, a and /3. Roll angle p is 
not decoupled, but will be treated as a discretely gridded input. This decoupling significantly promotes 
computational efficiency, which is crucial for on-line applications. The Hamiltonian function becomes: 


H (p, x, u) = yT ~Pi^~ v2c D a 2 a2 + (P 2 C L a COS <p - piV C D J a -p 2 ^VC Yl3 sin <p/3 (20) 

where p are the co-states of the value function: pi = an q p 2 = . This Hamiltonian is linear in thrust 

T and sideslip angle /3 and quadratic in angle of attack a. This structure allows for an efficient optimization 
routine over the inputs. The allowed ranges of the virtual inputs are: T m i n = 20546A/", T max = 41092IW, 
«min = 0°, a max = 14 -5° (no stall), 0 min /max = ±60°, (3 m[n / max = ±5°. The maximizers T, a and /? depend 
on the sign of the costates pi and p 2 - Recall that V > 0 , >0, Cd 2 > 0, >0, Cy p < 0. Due to the 

underlying physics, no sign changes for these parameters are to be expected in case of structural changes. 
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Define ft = p 2 a 2 lVClDa and a = amin ^ Qmax . Then the optimizing control inputs can be defined, 


according to the purpose. 

For invariance, the infmin optimization is achieved 
by means of the following inputs: 

• If pi > 0 then T = T min and 

if ft > a then a = a m [ n 

if ft = a then a = a m [ n or <a max 

if ft < a then a = <a max 

• If Pi = 0 then T e [Train! T max ] and 

if P 2 > 0 then a = <a m i n 

if P2 = 0 then QL G [c^min, ^max] 

if p 2 < 0 then a = <a max 

• If pi < 0 then T = T max and 

if p < G then a — ttm.in 
if Chmin ^ ft fS ^max then OL — ft 
if ft > a max then a = a max 


Alternative inputs are defined for the supmin opti- 
mization related to viability: 

• If pi < 0 then T = T max and 

if ft > ol then a = a m [ n 

if ft = a then a = <a m i n or <a max 

if ft < a then a = <a m ax 

• If Pi = 0 then f e [Train! T max ] and 

if P 2 < 0 then a = <a m in 

if P 2 = 0 then a G [a min ; a max ] 

if P 2 > 0 then a = a max 

• If pi > 0 then T = T m i n and 

if p ^ cr m i n then a — a^iji 
if Q^min A P A ^max then Qi = p 
if ft > Ol max then a = a max 


Besides, for the sideslip angle f3 holds: Besides, for the sideslip angle /3 holds: 

• if p 2 sin p > 0 then (3 = (3^^ • if p 2 sin p < 0 then j3 = /3 m i n 

• if p 2 sin ip = 0 then G [/ 3 min ; /3 max ] • if p 2 sin p = 0 then /? G [/3 min ; /3 max ] 

• if p 2 sirup < 0 then /3 = /3 max • if p 2 sin p > 0 then j3 = /3 max 

For the purpose of reachability analysis, the principle of duality can be applied for the infmin optimization 
inputs, as highlighted in Eq. 0. 


V. Results 

The first set of results focuses on motion in the plane of symmetry, assuming bank angle p = 0. Thereafter, 
motion leaving the plane of symmetry will be considered. 

Fig. [5] illustrates the definition of invariance, viability and reachability in the context of the RCAM 
example. All these properties are considered in the backward sense, meaning timesteps are taken backward 
in time. As a reference, the safe trim envelope K is defined as the area between the boundaries Tdnin = 
60m/s, Tmax = 100 m/s, 7 m i n = —10° and y m ax = +10° and is described by the implicit function : l (x) = 
min {xi — Emin, E max — aq, x 2 — 7min, 7max — X 2 } such that l (x) > 0 for x G K and l (x) < 0 for x ^ K. In 
Fig. |5(a)| the green area shows the invariant set X (T = 2s, K) as defined in Sec. |TTJ Monte Carlo analyses 
confirm the accuracy of the result. For initial conditions xo inside the invariant set Z(T = 2s, if), the 
whole state trajectory (ft (r, £, x, u) till the end point at T = 2s will remain inside the reference envelope if, 
independent from which inputs u (•) are applied to the aircraft. For initial conditions xo outside the invariant 
set X (T = 2s, if), there exists at least one input u (•) which will push part of the state trajectory (ft (r, £, x, u) 
towards the end point at T = 2s outside the trim envelope if. This complies with the aforementioned 
definition of invariant sets in Section El 

As can be seen in Fig. |5(b)| the viability set V (T = 2s, if) is somewhat less restrictive. Monte Carlo 
analyses confirm that, for initial conditions xq inside the viability set V (T = 2s, if), it is always possible to 
find at least one admissible input u (•) within saturation bounds which will keep the whole state trajectory 
(ft (r, £, x, u) till the end point at T = 2s inside the reference envelope if. On the other hand, outside the 
viability set V (T = 2s, if), the state trajectory (ft (r, £, x, u) is guaranteed to leave the reference envelope if 
within T = 2s, independent from which input u is applied. 

Finally, the concept of backwards reachability 7 Z(T = 2s, if ) is illustrated in Fig. |5(c)| For initial 
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conditions xo within the backwards reachable set 7 Z(T = 2 s,K), it is always possible to find at least one 
admissible input u(-) which will bring part of the state trajectory </> (r, £,x, u) towards the end point at 
T = 2s inside the trim envelope K. On the other hand, from outside the backwards reachability set 
7 Z(T = 2s, K), it is impossible for the state trajectory (p (r, £, x, u) to reach the reference envelope K within 
T = 2s, independent from which input u is applied. 

Fig. [6] shows the backwards evolution of invariance and viability sets over different discrete time steps 
of a total backwards time horizon T = 7s. At t = 7s, both sets are identical to the reference envelope K. 
Over the next few seconds, both sets start shrinking. While the invariant set keeps shrinking and eventually 
disappears altogether around t = 4.6s, the viability set converges to a steady state situation, which is reached 
around t = 3s. The difference between both can be analysed from their definition. For a time horizon long 
enough (in this specific case 4.6s), it is always possible to find at least one input within saturation limits 
which will push the state trajectory starting at whatever initial condition outside the reference envelope. The 
converging viability set towards a steady state shape indicates that there is a critical time instant (in this 
specific case 3s) for which at least one input can be found which keeps the state trajectories from the viable 
initial conditions within the reference envelope. Prolonging the time horizon won’t shrink this viability set 
further. This observation allows to make a time independent statement: the viability set for t = 3s and less 
is the set of initial conditions for which at least one admissible input can be found to keep the state trajectory 
within the reference envelope, irrespective of the length of the time horizon. Outside this viability set, the 
state trajectory is guaranteed to leave the reference envelope eventually at some point, independent from 
which input is applied. This is in contrast with the definition of the stable trim envelop^. As a consequence, 
it can be stated that under any circumstance the stable trim envelope should be situated inside the viability 
set. 

The backwards invariance, viability and reachability sets for the RCAM model over a time horizon of 
T = 2s are compared in Fig. [71 This set comparison confirms that invariance is a subset of viability, which 
in its turn is a subset of reachability, as was already claimed in Eq. ©• 

So far, only backwards integration has been considered. Besides, forward integration is important in 
the interpretation of the safe maneuvering envelope as elaborated in Section [U The only difference between 
both is the direction of the time horizon and the time steps: backwards reachability considers rearward 
timesteps, as has been discussed so far in this section, while forwards reachability involves forward time 
steps. Fig. [8] compares forward and backward reachability for the RCAM model over a time horizon T = 2s. 
Backwards reachability is the lower set and is consistent with previous figures, forwards reachability is the 
upper envelope. According to the definition of the safe maneuvering envelope provided in Section [I| the cross 
section between both in the center is designated as the safe maneuvering envelope. All flight conditions in 
this set can be reached from the stable trim envelope, and the stable trim envelope can be reached from this 
set as well. 

Despite the transparent definition of the safe maneuvering envelope provided previously, this set of flight 
conditions is somewhat too restrictive in the context of fault tolerant control and upset recovery. For this 
purpose, the backwards reachable set can be defined as the survivable envelope. This envelope corresponds 
to the set of flight conditions from which a safe recovery to a stable trim flight condition is still possible 
within the specified time horizon T after a disturbing failure, wind gust or other upset. As an illustration of 
this concept, the safe and survivable envelopes for the RCAM model over a time horizon T = 2s are given in 
Fig. [9] (dark and light green respectively), together with the reference envelope K which is marked by the red 
contour. The pilot can safely steer the aircraft arbitrarily as desired inside the dark green safe maneuvering 
envelope within 2s. In case of a failure, wind gust or other upset, it is possible to recover towards a safe 
equilibrium condition within 2s as long as the initial condition immediately after perturbation is situated 
inside the light green survivable envelope. However, given the admissible inputs between saturation limits, 
it is impossible for the pilot to steer the plane actively towards the light green area. In the next examples, 
focus will be placed on the survivable maneuvering envelope. 

Some damage scenarios are considered in Fig. [lOj and their influence on the survivable maneuvering 

a The trim envelope is the set of flight conditions for which an equilibrium can be reached, meaning that an admissible input 
within saturation limits can be found for these flight conditions which renders all values of the first order time derivative of the 
state vector equal to zero. Stability can be defined as static and dynamic stability. An equilibrium point is statically stable 
when a deviation from equilibrium, caused by a disturbance, will automatically generate a force or moment in the opposite 
direction with respect to the deviation. A dynamically stable trim point is a statically stable equilibrium point where the force 
or moment, automatically generated by the disturbance, is large enough to overcome and eliminate the disturbance. In this 
paper, a stable trim envelope is de facto the collection of flight conditions in a dynamically stable equilibrium. 
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A safe set at time = 2s 


An unsafe set at time = 2s 




(a) Backwards invariance with Monte Carlo Analysis 




(b) Backwards viability with Monte Carlo Analysis 




(c) Backwards reachability with Monte Carlo Analysis 


Figure 5. Backwards invariance, viability and reachability for RCAM model for time horizon T = 2s, including 
Monte Carlo Analysis 
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Time = 7s 


Time = 6s 


Time = 5s 


Figure 6. 


Figure 7. 

T = 2s 











Comparison of backwards invariance and viability for RCAM model over different time instants 


Comparison of reachability, viability and invariance 



Comparison of backwards invariance, viability and reachability for RCAM model for time horizon 
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Comparison of forwards and backwards reachability 



Figure 8. Comparison of forward (upper area) and backward (lower area) reachability for RCAM model for 
time horizon T = 2s 


Analsysis of safe set 



Figure 9. Safe (dark green) and survivable (light green) maneuvering envelopes for RCAM model for time 
horizon T = 2s. 
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envelope is investigated. Fig. |10(a)| considers a generic damage scenario where the total lift force is reduced 
by 20% and the total drag force is increased by 20%. As can be seen in Fig. |10(a)| the maneuvering envelope 
for the damaged aircraft shifts upward and expands slightly to the right. The physical explanation for this 
is the fact that the loss in lift force results in reduced climb capability but improved diving capability. As 
a consequence, the reference envelope can be reached from smaller negative, but larger positive flight path 
angles 7 . Moreover, increased drag leads to a faster deceleration capability, which makes the trim envelope 
reachable from slightly larger airspeeds within 2 s. 

As second damage scenario, an input restriction is defined where engine damage results in a 30% lower 
limit of the maximum thrust value. In Fig. |10(b)| it can be seen that only the left boundary of the envelope 
shrinks somewhat. This corresponds to the reduced acceleration capability due to the lower T max . The 
reference envelope can only be reached from slightly higher speeds in the lower range between 50m /s and 
60 m/s, compared with the nominal situation. 

Figure |10(c)| shows the cumulative effect of both damage scenarios, where both previous changes and 
shifts have been combined in the new survivable maneuvering envelope for the RCAM model over a time 
horizon of T = 2s. 




(a) 20% lift loss and 20% drag increase 


(b) 30% decrease of T max 



(c) 20% lift loss and 20% drag increase as well as 30% decrease 
of T'max 

Figure 10. Different damage scenarios for RCAM model for time horizon T = 2s 

So far, only motion in the plane of symmetry has been considered. By assuming that bank angle ip % 0, 
an additional dimension is added to the maneuvering envelope. The computationally most efficient way to 
do this, is defining a discrete set of bank angle values cp £ [0°;60°], and calculating the V, 7 maneuvering 
envelope for every value of cp individually. Due to symmetrical properties of the aircraft, only positive bank 
angles need to be considered. In case of asymmetric damage, the bank angle range needs to be adapted 
accordingly. The resulting 3D maneuvering envelope is given in Fig. [TlJ A more extensive analysis of these 
results will be given in section IVIi together with the incorporation of some robustness characteristics in the 
determination of the survivable maneuvering envelope. 
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Figure 11. Result of 2D optimization procedure over different bank angles </>, time horizon T = 2s 


VI. Robustness consideration 


In the optimal control framework elaborated previously, it is possible to incorporate robustness for un- 
certainties in the estimated values for the aerodynamic derivatives. These uncertainties are based on the 
standard deviations for the identified aerodynamic parameters: 

Cl. = Cl. + s L . • A L . with | A L# | ^ 1, (21) 

Cd. = Cd. + sd. • Ad. with |A£> # |^1, (22) 

Cy. = Cy. + sy. * A y. with |Ay # | ^ 1 (23) 

Based upon the required degree of conservatism, the standard deviations can be multiplied by 2 or 3 in 
order to increase coverage of the Gaussian distribution of the uncertain parameter (1 s 0# = 68%, 2s m 9 = 95%, 
3s.. =99%). 

The crux is to include the A’s as disturbances in the optimization function, they oppose the optimization 
over u. As a consequence, the invariance optimization formulation becomes!^ 


V 2 (x, t) = inf sup min l (0 (r, £, x, u (•) , A (•)))> 0 

u (')eW[t,T] A(-)e[-i,i] rG [ t ’ T ] 


(24) 


Similarly as previously, reformulation into an HJB PDE yields: 


dV? dV? 

— — (x, t) + min < inf sup — — (x, t) f (x, u, A) > = 0 

dt re[t,T] [ u(.)€M [t ,T] A(.)e[-U] dx I 


(25) 


where V 2 (x, T) = l (x) holds for backward integration and Vo (x, t) = l (x) applies to forward integration. 

For this purpose the Hamiltonian from Eq. flu can be rewritten, this time including parts independent 
of the inputs T, a or /3 but with some aerodynamic derivative(s): 


H (p,x,u) 


-Pitt-V 2 Cd 0 + P2^VC Lo costp - px^V 2 (C Da a + C D 2 a 2 ) +P 2 V^-C La cos pa + 
2m 2m 2m v " y 2m 

~P 2 V C Yfj sin tpP (26) 

2m p 
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It can be observed that the aerodynamic derivatives all appear linearly in an uncoupled way, which allows a 
similar procedure to solve the optimization as previously. By rewriting the Hamiltonian as a summation of 
terms, where each term is a multiplication of a variable involving a costate, a constant factor and a derivative, 
one can determine the sign of this factor, which consists of the predefined physical parameters: 


H( P,x,u) = 


-Pi o + p 2 guco Stfi C Lo - P1 + P2 ^Vaco Sl pC La + 


2 m 
>o 


pS_, 

2 m 


2 m 


>o 


>o 


>o 


“Pi ^~ v2(X C D a ~ P 2 sin ip/3~V C Y p 


ps_ 

2m 

>o 


(27) 


where it should be noted that p G [— 60°;60°] , a G [0°; 14.5°], /? G [—5°; 5°]. Furthermore airspeed V > 0 
and for the aerodynamic derivatives, it is known that Cd a > 0 , (7^ 2 > 0 , > 0, C Yf3 < 0. Due to 

the underlying physics, no sign changes for these parameters are to be expected in case of uncertainty or 
structural changes. 

Based upon this formulation, optimal control inputs for the aerodynamic derivatives can be defined as 
given in Table Q] where C. = C. max and C m = C. min . 


Table 1. Optimal control inputs for robustness against uncertain aerodynamic derivatives 


sign of costate 
Pi > 0 
Pi < 0 
P2 > 0 
P2 < 0 

P 2 sin p • P > 0 
P 2 sin p • (3 <0 


minimizer 

Cd 0 G Do •> G D a G ]j a 5 C D a 2 

Cd 0 — G-Do J GD a = C D(x , Co a 2 = Gd q 2 

C Lo =G Lo ,C La =G La 

C Lo =C Lo ,CL a =C L(x 
Or? = C Y p 

Cy p =C Yb 


maximizer 

Cd 0 — Gd () 5 = Gd cx 5 CD a 2 = Gd cm 2 

G Dq C Dq , C D a C D a 5 G G Y ) ^2 

C Lo =C Lo ,C La =C La 

c Lo = c Lo ,c La = c La 
= C Yp 


With this information, it is possible to create an entire “uncertainty band” around the envelope, however, 
here focus will be placed on the “worst-case” minimal size envelope. 

Fig. [12] compares the 3D envelopes with and without uncertainty, where two levels of uncertainty have 
been considered here, namely 10% and 20% of uncertainty on all aerodynamic derivatives. For the purpose of 
this example, identical ratios of standard deviations over nominal values have been defined for all derivatives, 
but the algorithm is capable to deal with individual standard deviations which can vary between the different 
aerodynamic derivatives. It can be clearly seen in Fig. [12] that larger degrees of uncertainty result in more 
significant shrinking of the envelope, since this is a “worst-case” minimal size envelope. 

Fig. [13] analyses the V, 7 maneuvering envelope for different values of bank angle p, including robustness 
for uncertainties of 10% and 20%. By comparing Fig. |13(a) , p~3(b) and |13(cj , it can be seen that larger bank 
angles have an influence on the climb capability of the aircraft. This is due to the physical principle that 
climb capability of lift force is provided through L cos p, which confirms a smaller decrease for smaller bank 
angles (up to p = 25° as shown in Fig. |13(b) ) but a much more significant change for larger bank angles as 
can be seen in Fig. |13(c) 


The final scenario to be considered, is the V, 7 maneuvering envelope for the RCAM model for different 
values of bank angle p, with and without generic damage scenario involving 20 % decrease in lift force and 
increase in drag force. The results are shown in Fig. [I4l and fl5l Fig. [Ml provides a comparison of the 3D 
envelopes without considering robustness. 

Fig. [15] includes robustness for 10% uncertainty and analyses the V, 7 maneuvering envelope for three 
different values of bank angle p, namely 0°, 25° and 60°. By comparing Fig. 15(a), |15(b)| and 15(c)| it 
can be seen that the damage scenario as well as larger bank angles both have a negative influence on the 
climb capability of the aircraft. This is again consistent with the underlying physical principles, as explained 
previously. For motion out of the plane of symmetry, climb capability is provided through L cos p. In this 
example, this capability is reduced by the non-zero bank angle p as well as by the reduced lift vector L, 
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Reachable region envelope of RCAM model at time = 0.00087s 



Y [deg] 


V [m/s] 


Figure 12. Comparison of 3D envelopes with and without uncertainty: nominal (green), 10% uncertainty 
(blue), 20 % uncertainty (yellow) 





(a) cp = 0° (b) ip = 25° (c) ip = 60° 


Figure 13. V, 7 maneuvering envelope of RCAM model for time horizon T = 2s for different bank angles and 
different uncertainty levels. Smaller envelope areas correspond to larger uncertainty bounds. 
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V[m/s] 


T [deg] 


Figure 14. Comparison of 3D envelopes with and without damage: nominal (blue), 20% change in lift and 
drag (green). No uncertainty included. 





(a) if = 0° (b ) if = 25° (c) if = 60° 

Figure 15. V, 7 maneuvering envelope of RCAM model without and with 20% change in lift and drag for 
time horizon T = 2s for different bank angles. 10% uncertainty included. Lower envelope areas correspond to 
nominal configuration. 
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resulting in an accumulated effect. This example is the combination of the situations in Fig. |10(a)| and [GO 
In case of a larger bank angle of cp = 60°, the decrease in climb capability due to loss of lift force is less than 
for zero bank angle, which is confirmed by comparing Fig. |15(a)[ |15(b) and 15(c)) 

Extensive Monte Carlo analyses have been performed in order to verify the accuracy of the boundary 
of the estimated survivable maneuvering envelope. These analyses have been based on the non-simplified 
aircraft model, ignoring the assumption that the aerodynamic angles a and {3 should be small. All these 
Monte Carlo analyses have confirmed that the results provided here are accurate and that the simplifications 
hold for the current ranges of the aerodynamic angles, namely a G [0°; 14.5°] (no stall) and /3 G [—5°; +5°]. 
This is an important conclusion which makes a relevant on-line safe maneuvering envelope estimation tool 
much more feasible. 


VII. Information required and output provided to other algorithms 

The research results presented here contribute to the VS ST project (Vehicle Safety Systems Technologies) 
and are planned to be integrated in the task package “Maneuverability Estimation & Upset Prevention”, 
involving the subtask activities “Dynamic Envelope Estimation & Protection (DEEP)” and “Constrained 
Trajectory Generation & Anticipatory Guidance (CTGAG)”. 

As a consequence, this maneuvering envelope estimation algorithm is interconnected with other algo- 
rithms, namely real time modeling, flight envelope protection, trajectory generation and anticipatory guid- 
ance. The first routine provides inputs to the algorithm, while its output is provided to the three latter ones. 
Figure [T6l illustrates these interconnections in a global overview. 



Figure 16. Overview of interconnection with other VSST components 

The algorithm for estimating the maneuvering envelope relies on some information which needs to be 
available. First, the aerodynamic derivatives and their standard deviations need to be estimated in real 
time, such that up to date information is available before and after failure. With this information, a stable 
trim envelope needs to be determined first. The values for derivatives and standard deviations, and the 
boundaries of the stable trim envelope serve as a-priori information for the algorithm discussed here. The 
output of the algorithm is then provided to the trajectory generation and anticipatory guidance algorithms. 
A more detailed overview of this can be found in Fig. [171 
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Tactical Flight Management System with Maneuverability Estimation for Upset Prevention 
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Figure 17. Detailed overview of task package maneuverability estimation and upset prevention 
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VIII. Conclusions and future work 


In this paper, a computationally efficient algorithm for estimating the safe maneuvering envelope of 
damaged aircraft has been discussed. The algorithm performs a robust reachability analysis through an 
optimal control formulation while making use of time scale separation and taking into account uncertainties 
in the aerodynamic derivatives. The safe maneuvering envelope is defined as the cross section between the 
forwards reachable and backwards reachable sets, which have been calculated starting from the stable trim 
envelope. Moreover, the backwards reachable set can be considered as the survivable maneuvering envelope, 
from where it is possible to bring the aircraft back to a safe trim condition after an upset due to damage, 
turbulence, a wake encounter etc. Damage scenario applications are included for two dimensional (in the 
plane of symmetry) and three dimensional examples. Results were found to be consistent with the underlying 
physical principles. This approach differs from others since it is physically inspired. This more transparent 
approach allows interpreting data in each step, and it is assumed that these physical models based upon 
flight dynamics theory will therefore facilitate certification for future real life applications. 

At this point, only the slow dynamics as specified in Fig. |2]have been considered. Future work will extend 
to the faster dynamics in higher bandwidths. It should be noted that only the low bandwidth dynamics have 
a full nonlinear input-output description, leading to a more complex scheme to determine the optimizing 
control inputs. The faster dynamics for both other ranges can be represented as a linear matrix input- 
output relationship, which simplifies the optimization scheme considerably. Moreover, the middle range is 
aircraft independent and does not need any robustness measures for uncertainties in aerodynamic derivatives. 
Another topic for further research is determining the bandwidth dependent time horizon values for every 
sub-part of the dynamics. 
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